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The experimental autoimmune encephalomyelitis (EAE) is an autoimmune disease of the central nervous 
system commonly used to study multiple sclerosis (MS). We combined clinical EAE phenotypes with 
genome-wide expression profiling in spleens from 150 backcross rats between susceptible DA and resistant 
PVG rat strains during the chronic EAE phase. This enabled correlation of transcripts with genotypes, other tran- 
scripts and clinical EAE phenotypes and implicated potential genetic causes and pathways in EAE. We detected 
2285 expression quantitative trait loci (eQTLs). Sixty out of 599 c/s-eQTLs overlapped well-known EAE QTLs and 
constitute positional candidate genes, including Ifiti (EaeT), Atg7{Eae20-22), Klrc3{eEae22) and Mfsd4{Eae17). 
A frans-eQTL that overlaps Eae23a regulated a large number of small RNAs and implicates a master regulator of 
transcription. We defined several disease-correlated networks enriched for pathways involved in cell-mediated 
immunity. They include C-type lectins, G protein coupled receptors, mitogen-activated protein kinases, trans- 
membrane proteins, suppressors of transcription {Jundp2 and Nrldl) and STAT transcription factors {Stat4) 
involved in interferon signaling. The most significant network was enriched forTcell functions, similarto genetic 
findings in MS, and revealed both established and novel gene interactions. Transcripts in the network have been 
associated with T cell proliferation and differentiation, the TCR signaling and regulation of regulatory T cells. 
A number of network genes and their family members have been associated with MS and/or other autoimmune 
diseases. Combining disease and genome-wide expression phenotypes provides a link between disease risk 
genes and distinct molecular pathways that are dysregulated during chronic autoimmune inflammation. 



INTRODUCTION 

by immune-mediated destruction of myelin sheaths and axons in 
A predisposition to develop a complex disease such as multiple the central nervous system, leading to progressive disability, 
sclerosis (MS) is regulated by numerous genetic variants that Despite recent substantial progress in deciphering genetic var- 
each contribute small effects (1). Clinically, MS is characterized iants that contribute to susceptibility (2), little is known about 
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the functional outcomes of these risk-associated variants, in part 
due to limitations in access of relevant human samples. 
In addition, identified risk alleles together explain only a fraction 
of disease heritability and variance (2). Additional risk variants 
conferring small effects may contribute to heritability of 
complex diseases, and the clustering of genes (below thresholds 
for significant association with MS) into functional networks 
supports this hypothesis (3). Unraveling the functions of suscep- 
tibility genes through identification of pathways enriched with 
risk genes can reveal mechanisms central in disease regulation. 

An animal model widely utilized to characterize the genetic 
basis and disease mechanisms of relevance for MS is experimen- 
tal autoimmune encephalomyelitis (EAE). Myelin oligodendro- 
cyte glycoprotein (MOG)-induced EAE in rats mimics many 
features of MS (4), including inflammation and demyelination, 
relapses and remissions and immune cell infiltration. Linkage 
analysis in experimental animal crosses can readily detect quan- 
titative trait loci (QTLs) related to clinical traits of complex dis- 
eases, and over 50 QTLs have been identified in EAE (5). Several 
genes underlying QTL effects in rats were positionally cloned 
and a number of them have been subsequently confirmed to regu- 
late human counterpart (6). However, it has been challenging to 
define single quantitative trait genes (7). 

Given the high heritability of variation in gene expression (8), 
identifying determinants of gene expression may give insights 
into pathogenic mechanisms of complex traits. The approach of 
mapping quantitative variation in gene expression was introduced 
in 2001 (9,10). This approach yields expression QTLs (cQTLs) 
(11), which influence expression of transcripts either in cis or in 
trans, where cw-acting eQTLs are located in close proximity of 
the target gene itself, while trans-acting eQTLs are located in a 
region distant from the gene it regulates. Technical artifacts 
excluded (e.g. hybridization differences) (12); c/s-regulatory 
effects can usually be mapped with high statistical significance 
and could be explained in most cases by a variation in DNA se- 
quence in the regulatory regions of the target gene (13). 
Cw-cQTLs that overlap trait QTLs constitute plausible candidate 
genes underlying the trait QTL effect. During the last decade, the 
use of genome-wide expression profiling combined with linkage 
analysis in segregating populations has identified genomic varia- 
tions that regulate complex traits in experimental models 
(6, 1 4, 1 5). In addition, the approach has been utilized to character- 
ize genetically driven networks of genes giving insights into path- 
ways and functions critical for the trait of interest (16). 

In this study, we combined genome-wide expression analysis in 
spleen from an experimental backcross (BC) between EAE- 
susceptible Dark Agouti (DA) and EAE-resistant Piebald Viral 
Glaxo (PVG) rat strains during the chronic phase of EAE with clin- 
ical EAE phenotypes and classical EAE QTLs. These two inbred 
rat strains have been extensively used in our laboratory to charac- 
terize EAE QTLs (5,17-27). We characterized several potential 
positional candidate genes for known EAE QTLs that provide a 
good base for further fianctional studies. Genome-wide expression 
analysis in the chronic stage EAE enabled correlation of transcripts 
not only with genotypes and to each-other but also with the clinical 
EAE phenotypes. We defined several disease-correlated gene net- 
works partially genetically regulated by loci that predispose for 
EAE. Some were enriched for pathways involved in cell-mediated 
immune mechanisms of relevance for EAE and MS, and also 
included genes or family members of genes associated with MS. 



RESULTS 

Overview of eQTLs in tlie clironic stage of EAE 

We used the eQTL approach to identify candidate genes and 
pathways that regulate EAE. This was achieved by combining 
clinical EAE phenotypes in a BC with whole-genome transcript 
expression analysis in splenic tissue from 150 BC male rats. 
Spleens were collected at day 35 after induction of EAE and ex- 
pression was measured using Affymetrix Rat Gene 1.0 ST 
Arrays. The evidence for linkage was tested between genotype 
and gene expression (27342 transcripts) to identify hereditary 
components and revealed a total of 2285 eQTLs with genome- 
wide significance of P < 0.05 (Table 1, Supplementary Mater- 
ial, Tables S 1 and S2). By introducing clinical traits as covariates 
that had been recorded during the clinical EAE experiments, we 
evaluated if the detection of eQTLs depends on disease inci- 
dence and severity. A majority of detected eQTLs did not 
depend on disease incidence or severity, as depicted in the 
Venn diagram (Fig. lA). Thus, we report only EAE status as a 
disease covariate in Table 1. An overview of cis- and trans- 
acting eQTLs across the genome is presented in Figure IB. 



Co-Iocalized trans eQTLs 

Trans-acting eQTLs (Fig. IB and Supplementary Material, 
Table S2) do not involve DNA variation in the expressed 



Table 1. Cis- and trans-eQTLs detected in EAE spleen day 35 post- 
immunization 



Chr 


Transcript 




EAE(a) 




EAE(i) 






cis 


trans 


cis 


trans 


cis 


trans 


1 


123 (20) 


264(12) 


106 (32) 


216(69) 


80 (26) 


71(15) 


2 


69 


196 


51 


142 


35 


46 


3 


23 


20 


19 


15 


17 


18 


4 


64 (47) 


8(6) 


62 (50) 


8(6) 


59 (50) 


13(7) 


5 


7 


17 


6 


13 


9 


20 


6 


45 


71 


39(8) 


57(4) 


37(9) 


105 (14) 


7 


14(2) 


5(1) 


13(1) 


3(1) 


8 


6(2) 


8 


25(9) 


12(4) 


27(2) 


24 


26(3) 


89 (3) 


9 


11 


4 


11 


6 


8 


8 


10 


21(17) 


24 (20) 


20(17) 


32 (29) 


14(11) 


22(21) 


11 


14 


35 


13 


45 


24 


509 


12 


10 


11(1) 


10 


10(1) 


10 


9 


13 


14(3) 


15(5) 


13(2) 


12(2) 


11(3) 


10(1) 


14 


13(3) 


17(2) 


10(3) 


9(3) 


9(3) 


6(2) 


15 


38 


49 


32 


31 


24 


8 


16 


31 


172 


25 


74 


25 


258 


17 


28(5) 


643 (141) 


24(7) 


408(102) 


10(2) 


158(31) 


18 


14 


20 


13 


17 


5 


7 


19 


15 


56 


13 


44 


11 


34 


20 


4 


34 


5 


38 


3 


14 


X 


16 


13 


17 


9 


15 


16 


Total 


599(106) 


1686(192) 


529(122) 


1213 (217) 


440(107) 


1427 (96) 



Number of cis- and /ra«.«-regulated transcripts for each chromosome and total 
(rows) for a selection of clinical phenotypes as covariates (column pairs). eQTLs 
were selected to have a logarithm of odds score > 2 (generated with the Haley- 
Knott regression model in R/qtl) and a genome-wide corrected P- value of < 0.05 
(generated with 1000 permutations). Numbers in parentheses refer to the subset 
of transcripts controlled by loci in EAE QTLs (see Supplementary Material, 
Table S5). Abbreviations: eQTL, expression quantitative trait locus; Chr, 
chromosome; Transcript, no covariate; EAE(a), incidence of EAE as covariate, 
additive model; EAE(i), incidence of EAE as covariate, interactive model. 
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EAE score day 35 
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Figure 1. Expression QTLs in chronic stage EAE. (A) Venn diagram indicating 
tlie number of overlapping transcripts wlien performing a whole genome scan for 
eQTLs using different clinical disease variables as covariates. (B) Scatter plot 
with the ,r-axis depicts the genomic position of the eQTL and the j'-axis depicts 
the genomic location of the target gene. Each circle represents a significant cis- 
or ^;'i7f!i-eQTL (no covariates used in the model). Blue lines represent EAE 
QTLs identified in the DA and PVG strain combination. The diagonal band indi- 
cates cis-eQTLs and the off-diagonal circles represent ?ra«s-eQTLs. Vertical 
fran^-eQTL bands reflect genomic regions that control many transcripts. 
eQTLs were generated with the Haley-Knott regression model in R/qtl; 
LOD > 2; genome-wide corrected P-value <0.05. 



transcript in question; instead, the transcript is regulated by other 
more distant genetic variations. Therefore, fraw^-regulation can 
denote a genomic location of a master regulator of transcription 
(16). Functional analysis using ingenuity pathways analysis 
(IP A) revealed a significant association of ?rawj-regulated genes 
with liver x receptor/retinoid x receptor activation, nuclear 
receptors involved in transcriptional regulation (28) and inactiva- 
tion of platelet activating factor (Supplementary Material, Tables 
S3 and S4). Vertical trans-bands of co-localized eQTLs reflect the 
genomic position of a regulator of multiple eQTL transcripts (29, 
30). Such bands were observed on rat chromosomes (RNO) 1, 2 



and 1 7, with the trans-band on RNO 1 7 overlapping an EAE regu- 
lating QTL, Eae23a, related to several clinical traits of complex 
disease (26). Although the trans-band on RNO 17 was not 
enriched for a functional pathway, we observed regulation of 
many small nucleolar RNA genes and spliceosomal small 
nuclear RNAs from this genomic region (Supplementary Mater- 
ial, Table S2). These have been reported to be involved in epigen- 
etic modifications (3 1) and the formation of spliceosomes (32), 
respectively. Additionally, we observed that genes in these fam- 
ilies negatively correlate with disease phenotypes. 



Candidate genes denoted by cis eQTLs 

A cM-acting eQTL directly identifies expression variations in a 
gene that in turn can regulate a physiological trait or molecular 
pathway. One of the most significant cw-eQTLs was the regula- 
tor of G-protein signaling 4 {Rgs4) (Table 2) on RNO 1 3 . We per- 
formed IPA analysis and recorded that Rgs4 and its most 
correlated transcripts associated with cell proliferation and 
migration, integrin and endocytosis signaling (Supplementary 
Material, Tables S3 and S4). Although Rgs4 does not reach a 
genome-wide significance in human GWAS, it is of interest 
that its family member Rgsl is associated with MS (33) and 
that the Rgsl4 gene region is associated with Crohn's disease, 
another disease with chronic inflammatory features (34). 

In the next stage, we combined eQTL analysis with classical 
EAE QTLs. Analysis of EAE QTLs that segregate between 
DA and PVG strains was performed in a larger set of 421 BC 
rats, encompassing 150 rats used for eQTL analysis, and all 
EAE QTLs are reported in Supplementary Material, Table S5 
(Stridh et al., unpublished data). A large proportion of the 
QTLs also showed evidence of linkage in a subset of 150 
animals used for eQTL analysis, albeit with less significance 
due to decreased power in a smaller sample. Therefore, we 
focused on all well-established EAE QTLs that have been iden- 
tified in more than one well-powered study (Table 3 , Supplemen- 
tary Material, Table S5), while emphasizing the QTLs identified 
in this BC (Table 3). 

We then examined the eQTL data for candidate genes under- 
lying EAE phenotypes. Of all eQTLs with a P- value <0.05, a 
total of 60 cw-eQTLs (Table 3) with LOD > 3.9 resided in pre- 
viously known EAE QTLs present in the same DA and PVG 
strain combination (Supplementary Material, Table S5), such 
as Ifitl (Eae7), Atg7 (Eae20-22), Klrc3 (eEae22) and Mfsd4 
(Eael 7), among others. These 60 cw-eQTLs were mainly asso- 
ciated with functional pathways of natural killer (NK) cell sig- 
naling and molecular functions of cellular growth and 
proliferation (Supplementary Material, Tables S3 and S4). A 
c/s-eQTL that maps to a disease QTL can be considered a 
likely causal gene underlying the disease QTL. 

Gene networks of inter-dependent genes that correlate 
witli clinical phenotypes 

Complex phenotypes are often the result of a response of mul- 
tiple functionally interacting genes. Combining genome-wide 
expression traits with clinical information enables the study of 
gene networks of inter-dependent genes that can be correlated 
with clinical phenotypes. With weighted correlation network 
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Table 2. Strong c/.«-eQTLs mapping outside of known EAE QTLs in DA/PVG. lavl crosses 



Chr 


LOD^ 


Mean 


Rat eQTL'' 


Genomic 


Gene symbol 


Gene name 




cis-eQTL 


Expression 


Probe set 


Location 








T7 T 
2.1 .1 


1 17.7 


1 mm A A z 


56/30912 


Fpr3 


Formyl peptide receptor 3 




12.0 


137.6 


1 A~7r\o inn 

10703706 


63694428 


Lilrb3 


Leukocyte immunoglobulin- like receptor, subfamily 


J 


13.2 


217.4 


1 AT ^ AO 

1071964K 


^ACA^OO/' 

79504886 


Zfp61 


Zinc finger protein 6 1 




lU.J 


3 ly.u 


iU /Uozy / 


yj /J4yj3 


Siglec5 


Sialic acid-binding Ig-like lectin 5 




12.9 


348.9 


1 A "71 -1 A"7T 

10724073 


139033461 


Art2 


t-Cell ecto-ADP-ribosyltransferase 1 




16.6 


78.5 


10724164 


1 CACAAAA/I 

i3y3yoyo4 


CnmalO 


Cholinergic receptor, nicotinic, alpha polypeptide 


} 


13.3 


59.8 


1 A "7 ^ ^ r\A A 

1071 1944 


1 nOAOT 1 CA 

iy8y67i30 


Ln'c27 


Leucine rich repeat containing 27 




1 1 "7 


ZJ i.O 


iU /2?syu4 


0 1 1 ooco^ c 
2 1 Jo28o38 


ivis4a / 


Membrane- spanning 4-domains, subfamily A, member 7 




lU. i 


1 C ^ 0 


1 mo on 1 0 
iU /2oy lo 


2 i41o2332 


Ms4a6a 


Membrane- spanning 4-domains, subfamily A, member 6 




1 0.3 


y /.y 


1 m 1 A 1 n/^ 


0 1 c cnn/i 1 
2i33yUoio 


Faml 11a 


Family with sequence similarity 111, member A (NP_001 102633. 1) 


2 


14.6 


413.1 


10s23903 


173683154 


Gucylb3 


/~> 1 J- 1 11 1_1 1 J- T 

uuanylate cyclase 1 , soluble, beta 3 


2 


28.6 


145.2 


1 AO ^ zi ^ A A 
lOo 16144 


\ 1 Z A TA") 1 A 

175479320 


Sfrp2 


Secreted frizzled-related protein 2 


2 


12.0 


263.3 


10s24357 


1809741 1 1 


Mstol 


Misato homolog 1 


2 


12.7 


1 14.7 


1 Ad A ^ ^ ^ 

1082461 1 


1 O"^ A A AAO 

182644098 


Slc27a3 


Long-chain fatty acid transport protein 3 


2 


10.7 


221.3 


1 AO 1 A/'AA 

10819690 


A A ^ '^'^ A 1 1 

24432341 1 


Mcoln2 


Mucolipin 2 


3 


10.2 


78.6 


1 AO A All ^ 

1084433 1 


1 1 C 1 1 /I AT 

1151 1403 


Lcn2 


Lipocalm 2 


■J 

D 




iDo.y 


1 no A mzn 
iU84 /y3 / 


nn 1 T/inco 
yul /4y32 


Prrg4 


Transmembrane gamma- carboxy glutamic acid protein 4 


■3 
J 




VLv.L 


iU84oo32 


1 nA^ /I Qo 1 ^ 
iUo348o 13 


j!>ptDn3 


Spectrin, beta, non- erythrocytic 5 


A 

4 


1 1 A 

1 i .4 


400.7 


iU8oo3U / 


i /i884j 1 / 


Art4 


Ecto-ADP-ribosyltransferase 4 


4 


10.0 


334.6 


1 AO cm T "7 

10859337 


174009860 


Pde6h 


Phosphodiesterase 6H, cGMP- specific, cone, gamma 


5 


14.9 


136.4 


1 AO"7C /lie 

10873425 


1 /I OAT zn") 

248y3373 


KCjU 1309085 


Similar to r23N ly.y 


c 

J 


1U.2 


1 1 /I T 

1 14. / 


iU8o / /oi 


1 0 1 CO 1 1 1 


Mmp 1 6 


Matrix metalloproteinase- 1 6 


c 

J 


inn 


/;i o 
oz.o 


1 noT? 1 1 n 
iU8 / / 1 JU 


/ /UiU /84 


Ptgrl 


Prostaglandin reductase 1 


c 
D 


lU.o 


350. / 


1 no/^m o 0 
iU8Dy2oa 


/ /y24233 


Snx30 


Sorting nexin family member 30 


6 


10.7 


1359.5 


1 AOAI/'d 

10892653 


142735827 


Igh-6 


Immunoglobulin heavy chain 6 


6 


25.7 


267.3 


1 AO A1 ^ <n 

108y2662 


1 A 1 AT7<1 OA 

142y77650 


IgG-2a 


uamma-2a immunoglobulin heavy chain 


7 


1 1.6 


260.5 


1 AAA A C ^ f\ 

10904539 


1 12913724 


Ly6k 


Lymphocyte antigen 6 complex, locus K 


7 


16.1 


536.7 


1 AAA A C A"7 

ioyo45y7 


1 ^ 1 Ai A A nn 

1 13434499 


Ly6a 


Lymphocyte antigen 6 complex, locus A (predicted) 


7 


15.0 


509.1 


1 AOAO 1 A/' 

10898196 


121816012 


Mppedl 


If J-„11 1 1 J- J 1- ■ ■ 1 

Metal lophosphoesterase domain containing 1 


o 

y 


16.3 


1 T7 0 


1 nno m o o 
iUy2y28o 


-rnnn i /i o c 
/yuy 1423 


Seipine2 


Serine (or cysteine) proteinase inhibitor, clade E 


1 1 


lU.D 


1A Z 


1 mc \ A1A 
iU /3 i4j4 


(-.(^1 Z.T1 1 /I 

ooi3 / / 14 


Csta 


Cystatin A (stefin A) 


1 1 


12.6 


76.6 


1 mc /I on /L 

10734876 


TAz: Z") AAT 

70633yy3 


LUL6543O0 


Similar to SMP3 mannosyltransferase 


12 


15.8 


109.3 


1 A "7 /■ 1 1 C >1 

10762254 


T/'OA/'OI A 

36896819 


Oaslk 


2 —5 oligoadenylate synthetase IK 


13 


29.0 


158.8 


10769672 


85533882 


Rgs4 


Regulator of G-protein signaling 4 


13 


10.2 


301.6 


10765469 


86096465 


Sh2dlbl 


Similar to EWS/FLIl activated transcript 2 


14 


24.2 


135.0 


10776325 


33403956 


LOC498350 


Similar to testicular haploid expressed gene product isoform 2 


14 


10.3 


45.5 


10778820 


109197351 


Ccdc85a 


Coiled-coil domain containing 85A 


16 


17.5 


229.3 


10787757 


22852487 


Csgalnactl 


Chondroitin sulfate W-acetylgalactosaminyltransferase 1 


16 


26.5 


409.1 


10789670 


84885522 


Lig4 


DNA ligase 4 


18 


15.4 


140.2 


10804396 


40914028 


Cdol 


Cysteine dioxygenase 1 , cytosolic 


X 


14.4 


236.7 


10936742 


21117225 


GPR34 


G-protein-coupled receptor GPR34 (predicted) 


X 


16.9 


73.6 


10937013 


26697034 


Kcndl 


Potassium voltage-gated channel, Shal-related subfamily, member 1 


X 


26.8 


106.0 


10932773 


35143324 


Chrdll 


Kohjirin 


X 


11.8 


68.2 


10939319 


122209367 


Armcx6 


Armadillo repeat containing, X-linked 6 


X 


12.8 


209.6 


10939498 


127688736 


LOC678934 


Similar to CG30327-PA 



eQTLs were selected to have a logarithm of odds (LOD) score > 10 and a genome-wide corrected P-value <0.05 (generated with 1000 peimutations). C«-eQTLs 
were considered those having no genetic marker between the peak of the linkage score and the chromosomal region coding the transcript. Genes in bold are 
differentially expressed between DA and PVG in day 7 post-immunization ex vivo LN cells and/or in MOG re-stimulated conditions (49). 
Abbreviations: QTL, quantitative trait locus; LOD, logarithm of the odds ratio; eQTL, expression quantitative trait locus; Chr, chromosome. 
°LOD score of cu-eQTL for probe sets generated with the Haley— BCnott regression model, no covariates, in R/qtl. 
''Average probe set expression level. 

°Probe ID annotation from Affymetrix Rat Gene LOST array. 
''Genomic location of the transcript, start, in mega base pairs. 



analysis (WGCNA) (35), we identified eight gene networks that 
significantly correlated with EAE phenotypes (including sus- 
ceptibility and severity phenotypes, weight loss and anti-MOG 
IgG titers). We performed IP A of transcripts for each separate 
gene network to identify their functional properties. Significant 
associations with canonical pathways were discerned for six of 
the eight identified gene networks (Supplementary Material, 
Table S6) and four of these are described below. 



A gene network enriched for T cell functions shows strong 
correlation with EAE 

The most striking gene network A (Fig. 2, Supplementary Mater- 
ial, Tables S6 and S7) associated with molecular functions 
including T cell-mediated immune mechanisms (Tables 4 
and 5), which have also been implicated in MS (2). This 
network gave the strongest positive correlation with EAE 
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Table 3. Os-eQTLs mapping to EAE QTLs in DA/PVG. lavl crosses 



E AE QTL 


Chr 


EAE" 


lod'' 


Mean" 


Rat** 


Genomic'' 


Gene symbol 


Gene name 






QTL 








eQTL 












Position 


c«-eQTL 


Expression 


Probe Set 


Location 






Eae29 


1 


0- 


25 


20.8 


34.5 


lu/Olozu 


5473 


RCjD15641 10_preaictea 


Similar to ptitative pheromone receptor 
(predicted) 


Eae29 








0 ^ 

O.J 




lU / IojZo 


Z4Moi 


HJL.3ojUZy 


Similar to vomeronasal 2, receptor, 1 


Eae29 








24.8 


80.6 


Iu701o6o 


1498164 


LUL67874U 


Similar to vomeronasal 2, receptor, 1 


Eae29 


1 






42.4 


34.5 


10701671 


1503580 


LOC687363 


Similar to vomeronasal 2, receptor, 1 


Eae29 








z4. J 


A 1 A 
41 .4 


lU / lOJOZ 


1 jZ4 /fsy 


T r^(~">0AO0A 

i^OL,zooyoo 


Putative pheromone receptor Go-VN13C 


Eae29 








15.8 


64.7 


lu71o56o 


^ A ~7no 

1524798 


T /~\ "> O z; A O z; 

LUL28oy8o 


Putative pheromone receptor Go-VN13C 


Eae29 








38.9 


63.6 


10701674 


1636457 


LOL308240 


T T J-1 J- ■ IT 1 A O 'T ^ A 

Hypothetical LOL308240 


Eae29 








20.5 


68.0 


10701684 


1 6943 1 5 


RCjD 1 565235_preaictea 


Similar to Retinoic acid early inducible protein 1 


Eae29 








14. U 


60.7 


1 f\n 1 A Am 
iU / iooUi 


2121666 


KLrU 1 JO jZ3 j_preaictea 


Similar to Retinoic acid early inducible protein 1 


EaeSO 








7.5 


A on A 


1 f\n 111 nn 

lu / 1 izyy 


1 Qni OA 10/1 

1 o /3yDl o4 


Itgax 


Integrin, alpha X 


EaeSO 








O.Z 


388.2 


1 f\n 1 1 A A/1 
lU / 1 1od4 


lyuyo /oj / 


Acadsb 


Short/branched chain specific acyl-CoA 
dehydrogenase 


EaeSO 


1 






5.7 


67.0 


10726371 


193384041 


Adam 1 2_predicted 


A disintegrin and metallopeptidase domain 12 


Hue/ 




238- 


258 


4.0 


AOT 1 
OZj. 1 


iU / i4VU / 


llOAHQ 1 AO 


Tfif 1 


Interferon-induced protein with tetratricopeptide 
repeats 1 


Eae24 


4 


59- 


75 


1 1.9 


68.2 


1 AO C /I /I A£ 

10s5440o 


£ 1 "7 1 /I 1 AA 

61714290 


A 1 . 1 1_ O 

AKrlb8 


Aldo-keto reductase family 1 , member B8 


Eae24 


4 






7.5 


200.2 


10861946 


C /' 1 A OAO 

65674808 


D630045J12R1K 


t T^'V TCI Tl "V T/^T^A A AA A A y1 A T A 1 1 

EN SRN 0 1 0000004039 1 


Eae24 


4 






5.6 


221.2 


10b546J7 


AATITJ/IO 

66232348 


Clec21 


c- 1 ype lectin domain iamily 2, member L 


Eae24 


4 






8.5 


11 .6 


10854942 


6901 1941 


RGD 1 560283_preaictea 


Similar to trypsinogen 8 (predicted) 


Eae24 


4 






13.4 


91.4 


10862184 


69022614 


RGD1560283_predicted 


Similar to trypsinogen 8 (predicted) 


Eae25 


4 


75- 


83 


3.9 


437.9 


10855512 


78343646 


Mpp6 


Membrane protein, palmitoylated 6 


hae2u-22 


4 


133- 


159 


6.7 


523.6 


1 AO C "7AC A 

10057950 


1 C AO 1 A "7AA 

150814706 


Atg7 


Autophagy-related 7 (yeast) 


Eae20-22 


4 






21.1 


89.4 


10865186 


156705371 


RGD1310710_predicted 


Similar to RIKEN cDNA 270009 1N06 
(predicted) 


17 -1/1 T» 

hae2u-22 


4 






8.5 


50.5 


1 AOAcn 1 

lOoojJzl 


1 Z OA A 1 O T C 

158y0l535 


RCjD 1560851 _preaictea 


Similar to mKIAA1238 protein (predicted) 


Eae20-22 


4 






10.9 


141.6 


10858581 


159745231 


Clec4b2 


c-Type lectin domain family 4, member B2 


eEae22 


4 


162- 


172 


6.9 


167.4 


1 AOCAAAA 

10859090 


1 /■/'AAAOAA 

166090800 


T /^/^/"OAOAA 

LOL689800 


Similar to osteoclast inhibitory lectin 


eEae22 


4 






9.3 


45.9 


1AOCA1 AC\ 

10859149 


166857707 


Klrel 


Killer cell lectin-like receptor, family E, member 


eEae22 


4 






5.5 


86.2 


10859162 


166857707 


Klrel 


T^'1l 111 I'l J- C 'IT^ 1_ 

Killer cell lectin-hke receptor, lamily E, member 


eEue22 


4 






3. J 


01/1 y 


1 nO^Q 1 A/1 

lUojy 1d4 


1 AAO ^ TTm 
lODOJ / /U / 


Ti^l^^l 1 

isJral 


Killer cell lectin-like receptor, subfamily D, 
member 1 


eE(ie22 


4 






J.O 


lU.o 


1 nO'^Q 1 AH 

lUojy lou 


1 AAO ^ TTm 
lOOOJ / /U / 


f l^a 1 


Killer cell lectin-like receptor, family E, member 


eEae22 


4 






4.9 


804.2 


10866041 


166914788 


Klrkl 


Killer cell lectin-like receptor subfamily K, 
member 1 


eEae22 


4 






13.2 


161.6 


IUsdoUjo 


1 A An C A A 1 O 

looyjool o 


Klrc2 


Killer cell lectin-like receptor subfamily C 


eEae22 


4 






17.6 




IUodoUjZ 


1 A An C A A 1 O 

looyjool o 


Klrc3 


Killer cell lectin-like receptor subfamily C, 
member 3 


eEae22 


4 






4.Z 




1 nCAAHA 1 


1 AAno 1 1 n 
looyo 133Z 


isJrcl 


Killer cell lectin-like receptor subfamily C 


eEae22 


4 






"7 /I 

/.4 


y /.u 


1 f\Q i^t.r^n c 
lUSDOU /o 


1 A'rnn 

lo /Uy 1333 


Klri2 


Killer cell lectin-like receptor family I member 2 


eEae22 


4 






1 1 ^ 


joU.z 


1 no AA 1 fl 1 

lUoOolUl 


1 ATIO^ 1 1 1 


KLrUi jo3i iu_preuictea 


Similar to immunoreceptor Ly49si3 


eEae22 


4 






/I 1 
4. 1 




1 no AA 1 1 A 


1 ATI ATI C n 

lo /3o /3 jy 


LL)l,oo4Ujy 


Similar to immunoreceptor Ly49sil 


eEae22 


4 






4.7 


229.4 


lOoDolzO 


1 £ "7 /I /I AO 1 A 

167449814 


Ly4ysil 


Immunoreceptor Ly49sil 


eEae22 


4 






10.4 


4UZ.O 


I nOAA 1 T3 
lUoDOlzi 


1 AT/l 0 m 1 o 
lo /453U1 o 


Klra22 


Killer cell lectin-like receptor subfamily A, 
member 22 


eEae22 


4 






1 1.4 


820.8 


1 AO AA 1 /IT 

1Uo6o142 


1 ATCAAIOA 

167500280 


T /"lO 1 AAT A /nc 1 

LUL10U364751 


Immunoreceptor Ly49si3-like 


eEae22 


4 






12.1 


472.3 


lAOAAl A A 

10bdo144 


1 ^1 Zr\ A CAT 

167504502 


T 1 AAT <1 /IT C 1 

LUL10U364751 


Immunoreceptor Ly49si3-like 


eEae22 


4 






10.7 


OlZ.J 


1 r\Q ^ Af\ 
lUoDol4U 


\ t^n zzic\Ai 
lo / JJ3U43 


T c\r^ \ (\r\i i-. An z \ 
HJL.1UU3o4/j1 


Immunoreceptor Ly49si3-like 


eEae22 


4 






10.7 


25.4 


1 AOA£ 1 AT 

10866163 


167y75702 


T ^,AC\„ A 

Ly4ys4 


Ly49 stimulatory receptor 4 


eEae22 


4 






19.6 


25.5 


10866186 


168230133 


Ly49i3 


Immunoreceptor Ly49i3 


eEae22 


4 






4.3 


172.4 


10866167 


168230133 


T ^.AC\„ A 

Ly4ys4 


Ly49 stimulatory receptor 4 


eEae22 


4 






6.5 


94.7 


10866215 


168575551 


NP_00 1009494.1 


Ly49 stimulatory receptor 7 


eEae22 


4 






1 A Q 
14. 0 


jZ.Z 


1 nOAATlA 

lUoOoZ JO 


1 A0"7/1 "7 1 1 n 

loo /4 / 1 lU 


isJraj 


Killer cell lectin-like receptor, subfamily A 


eEae22 


4 






7.2 


58.7 


10866243 


168794389 


Ly49i7 


Immunoreceptor Ly49i7 


eEae22 


4 






20.6 


59.8 


10866359 


170271493 


LOC689869 


Similar to Taste receptor type 2 member 140 


Eae_mo4 


4 


175- 


187 


9.1 


55.5 


10867020 


186999579 


RGDl 56085 l_predicted 


Similar to mKIAA1238 protein (predicted) 


Eae3 


10 


15- 


45 


10.8 


423.5 


10742431 


35933369 


Rufyl 


RUN and FYVE domain containing 1 


EaelSa 


10 


55- 


67 


4.7 


119.1 


10734740 


55270065 


Pilc3r6 


Phosphoinositide-3 -kinase, regulatory subunit 6 


EaelSa 


10 






4.3 


75.7 


10744254 


56692319 


Tnkl 


Non-receptor tyro sine -protein kinase TNKl 


Eael2 


10 


85- 


108 


8.2 


189.0 


10738451 


90602484 


Tmeml06a 


Transmembrane protein 106 A 


Eael2 


10 






4.9 


5793.3 


10774359 


104545644 


Rpl38 


Ribosomal protein L38 


Eael2 


10 






15.9 


238.8 


10748857 


105084291 


RGD1561778 


Similar to dendritic cell-derived 
immunoglobulin(Ig)-like receptor 1 


Eaell 


13 


38- 


55 


6.2 


138.9 


10763751 


43752437 


Fcamr 


Fc receptor, IgA, IgM, high affinity 



Continued 
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Table 3. Continued 



EAEQTL Chr EAE" LOD'' Mean" Rat'' Genomic" Gene symbol Gene name 

QTL eQTL 
Position c«-eQTL Expression Probe Set Location 



Eael7 


13 






4.3 


270.9 


10767565 


44910207 


Mfsd4 


EaelO 


14 


0- 


-20 


22.7 


236.9 


10771171 


5807837 


LOC683128 


EaelO 


14 






7.4 


142.3 


10771190 


5965433 


Abcg313 


Eae23b 


17 


57- 


-66 


12.7 


69.2 


10795634 


63355581 


4921524L21Rik 


EaelS 


18 


69- 


-87 


6.0 


116.7 


10805571 


81524266 


Timm2 1 



Major facilitator superfamily domain-containing 
protein 4 

Similar to guanylate binding protein family 
ATP-binding cassette, sub-family G, member 3 
family member 

Translocase of iimer mitochondrial membrane 
21 homolog 



The first column gives the name of the previously mapped EAE QTL (See Supplementary Material, Table S5 for related publications). The QTLs given in bold 
displayed linkage in 42 1 BC rats, encompassing 150 rats used for eQTL mapping. eQTLs were selected to have a logarithm of odds (LOD) score > 3 .9 and a 
genome-wide corrected /"-value < 0.05 (generated with 1 000 permutations). 0.s-eQTLs were considered those having no genetic marker between the peak of the 
linkage score and the chromosomal region coding the transcript. Genes in bold are differentially expressed between DA and PVG in day 7 post-immunization ex vivo 
LN cells and/or in MOG re-stimulated conditions (49). 

Abbreviations: eQTL, expression quantitative trait locus; QTL, quantitative trait locus; Chr, chromosome; LOD, logarithm of the odds ratio. 
^Position of EAE QTL in mega base pairs. 

''LOD score of c«-eQTL for probe sets generated with the Haley-Knott regression model, no covariates, in R/qtl. 
"Average probe set expression level. 

''Probe ID aimotation from Affymetrix Rat Gene 1 .0 ST array. 
"Genomic location of the transcript, start, in mega base pairs. 



susceptibility and severity (Pearson correlations: EAE p = 0.23, 
p = 0.004; MAX p = 0.23, P = 0.006; SUM p = 0.32, P = 
0.00006; s35 p = 0.26, P = 0.001; ONS p = -0.21, P = 0.01; 
DUR p = 0.3 1 , P = 0.000 1 ; WL p = 0.32, P = 0.00007). Tran- 
scripts in gene network A (Supplementary Material, Table S7) 
included multiple genes associated with autoimmune diseases 
iCd6, Etsl, Tcf7 and Themis) (2, 36-38). The Cd6MS suscepti- 
bility allele is associated with alterations in T cell proliferation 
(39) and the transcription factor Etsl participates in important 
aspects of early thymocyte development (40, 41). Several 
genes, including Lefl, Lck, Crtam and Itk, have previously 
been linked to functions of the adaptive immune system. In the 
network, the cw-regulated major facilitator superfamily (MPS) 
domain containing protein 4 (Mfsd4), overlapping an 
EAE-regulating QTL on RNO 1 3 that was initially identified in 
a (LEW X PVG) F2 intercross (42) had a strong (-0.05) 
albeit peripheral membership. Lefl, co-regulated by two loci 
on RNOl and RNO 18, was the strongest member of this gene 
network ( — 0.21). Several members of the network are regulated 
from EAE QTLs, Art2 (EaeSl), Crtam {Eael 7 and Eae23a) and 
ResplS {Eae23a), showing a genetic regulation of the network 
by loci that predispose for EAE. 

Gene network associated with G protein coupled 
receptor signaling 

Gene network B (Fig. 3, Supplementary Material, Tables S6 and 
S8) con^elated negatively with severity phenotypes (MAX 
p=-0.19, P=0.02; SUM p=-0.24, /'= 0.003; s35 
p=-0.19, P=0.02; DUR p=-0.16, P = 0.05; WL 
p= —0.27, 0.0005). Canonical pathways associated to 
the gene network included LPS/IL-1 -mediated inhibition of 
RXR function and G protein coupled receptor signaling (Supple- 
mentary Material, Table S6). Several transcripts in the network 
are family members of genes associated with MS (Supplemen- 
tary Material, Table S8). These include C-type lectins, G 



protein coupled receptors, mitogen-activated protein kinases, 
transmembrane proteins and solute carriers (SLCs) which are 
proteins involved in both inflammation and immunity, as well 
as in signal transduction that leads to cellular responses. 
Jundpl and Nrldl are in regions associated with Crohn's 
disease (34) and are both involved in transcriptional regulation 
(43, 44). The gene network is partially genetically regulated by 
Eael 7, Eae23a, Eae29 and EaeSO. 

Genetic regulation of gene network associated with IFN 
signaling pathways 

Gene network C (Fig. 4, Supplementary Material, Tables S6 and 
S9) con-elated with duration of EAE (p = 0.18, P = 0.02) and 
indicated a genetic regulation by loci that predispose for EAE 
{Eae3, Eae7, EaelSa, EaelSb and Eae20-22). IPA revealed a 
high significance for IFN signaling pathways (Supplementary 
Material, Table S6) that are important in neuroinflammation. 
Network C also included genes and family members of genes 
associated with MS or other autoimmune diseases; zinc 
fingers, STAT transcription factors {Stat2 and Stat4) (45), 
genes involved in the innate immune response (Irgm) (46) and 
interferon regulatory factors {Irfl, Irf7 and Ir/9) (45, 47). 
IFN-induced downstream signaling molecules, several of 
which associated with MS (48), are evident in the network 
such as proteins induced by interferons (Oas), interferon- 
induced guanylate-binding proteins (Gbp) and a number of 
interferon-induced proteins. 

Furthermore, network D (Supplementary Material, Table S6) 
correlated with serum anti-MOG total IgG levels day 35 p.i. (p = 
0.24, P = 0.003) and was functionally enriched for the regula- 
tion of p70S6K, a serine/threonine kinase in the PI3 kinase 
pathway. 

In conclusion, we identified several gene networks involved in 
specific immune-related pathways that are presumably of im- 
portance or crucial for both EAE and MS pathogenesis and 
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Figure 2. Gene network enriched for T cell functions shows strong correlation with EAE. Expression traits with neighboring locations are drawn next to each other. 
eQTLs are indicated by their rat chromosome location on a second line and c«-eQTLs are indicated in white. eQTLs within EAE QTLs are indicated by a grey box with 
the QTL name listed in parentheses. Predicted hub genes are presented in yellow and the black arrow heads indicate a directed interaction. Plus and minus indicate 
association to EAE phenotype. Genes with a higher disease association are placed to the right, whereas causative genes are placed toward the left. 



which indicate partial genetic regulation by loci that predispose 
for EAE. Additionally, the co-expression analyses provided 
novel insights in the cause-effect interactions of several genes 
that have not previously been linked to immune related path- 
ways. Furthermore, the approach of integrating clinical data in 
the analyses distinguished molecular responses important at dif- 
ferent stages of disease. For instance, we defined networks of 
interacting genes that correlated with duration and IgG levels, re- 
spectively. Any genetic difference in individuals that could 
modulate the defined molecular pathways could potentially in- 
fluence susceptibility to EAE and MS. 



DISCUSSION 

The eQTL approach can facilitate identification of candidate 
genes (6,14,15) and in the current study our whole-genome ap- 
proach identified genetic differences between two inbred rat 
strains that contribute to the expression of individual genes. 
We detected a number of c/s-eQTLs that overlap established 
EAE QTLs that warrant further biological characterization. 



The co-localization of a cw-eQTL and an EAE QTL makes it 
an attractive novel positional candidate for the causal relation- 
ship between the disease phenotype and the regulation of the 
transcript. We identified several cw-regulated positional candi- 
date genes in EAE QTLs including Ifitl (Eae7), Atg7 
(Eae20-22), Klrc3 (eEae22) and Mfsd4 (Eael 7). We have previ- 
ously described these genes as being differentially expressed 
between DA and PVG in ex vivo LN cells and/or in MOG 
re-stimulated lymphocytes at day 7 p.i. of EAE (49). Ifitl is 
induced in response to stimuli such as LPS, IL-1 and TNF 
(50,51) and has been described among genes that could best 
predict the response to IFN-|3 treatment in MS patients (52). 
Atg7 is involved in macroautophagy, an intracellular pathway 
that regulates quantity and quality of organelles and proteins. 
Autophagy has recently been linked to both innate and adaptive 
immunity (53). Klrc3, encodes an activating NK cell receptor, 
and NK cells are suggested to confer a disease -protective role 
in MS and EAE (54-56). Mfsd4 is a transmembrane transporter 
of the MFS. MFS proteins are the largest group of secondary 
carriers in the cell and transport small solutes by using chemios- 
motic ion gradients (57,58). The availability of whole-genome 
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Table 4. Top canonical pathways associated with Network A 



Ingenuity top canonical pathway 



Ratio" Genes'' 



f-value" 



Corrected 
F-value'' 



T Cell receptor signaling 
iCOS-iCOSL signaling in T helper cells 
CD28 signaling in T helper cells 
Calcium-induced T lymphocyte 

apoptosis 
Role of NFAT in regulation of the 

immune response 



1 1/109 CD247,CD3G,CD28,LCK,CAMK4,PRKCQ,CD3E,CD8A,CD3D,CD8B,1TK 

10/122 CD247,CD3G,CD28,LCK,CD40LG,CAMK4,PRKCQ,CD3E,CD3D,1TK 

9/132 CD247,CD3G,CD28,LCK,CAMK4,PRKCQ,CD3E,CD3D,1TK 

7/70 CD247,CD3G,LCK,CAMK4,PRKCQ,CD3E,CD3D 

9/200 CD247,CD3G,CD28,LCK,CAMK4,PRKCQ,CD3E,CD3D,1TK 



2.1E- 
1.3E- 
1.6E- 
1.3E- 



15 
13 
11 
10 



3.1E-13 
9.6E-12 
7.9E-10 
4.8E-09 



4.8E-10 1.4E-08 



Transcripts in Network A were uploaded into the ingenuity pathways application (Ingenuity Systems, www.ingenuity.com). The software was used to identify the 
most significant pathways. 

"Ratio refers to the number of molecules in a given pathway that meet cutoff criteria, divided by the total number of genes that map to the canonical pathway. 
''Genes in a given pathway that meet cutoff criteria. 

"Significance was determined from a single test /"-value calculated using the right-tailed Fisher's exact test. 

''Adjusted significance was determined from multiple test-corrected P-values using the Benjamini-Hochberg coiTection (top five categories reaching a coiTected 
statistical significance of <0.05 are shown). 



Table 5. Molecular functions associated with Network A 



Cellular and molecular functions 


Molecules" 


Representative subgroup 


Molecules'' 


f-value" 






P-value subgroup'' 


Cellular function and maintenance 


25 


Lymphocyte homeostasis 


23 


1.3E-23- 


1.2E 


-02 


1.3E-23 


Cellular development 


27 


T cell development 


22 


3.0E-23- 


1.2E 


-02 


5.9E-23 


Cell-to-cell signaling and interaction 


23 


Activation of T lymphocytes 


15 


3.6E-16- 


1.3E 


-02 


3.6E-16 


Cellular growth and proliferation 


26 


Proliferation of lymphocytes 


18 


3.8E-15- 


1.2E 


-02 


9.4E-15 


Cellular movement 


18 


Lymphocyte migration 


13 


5.0E-13- 


1.2E 


-02 


5.0E-13 



Transcripts in Network A were uploaded into the ingenuity pathways application (Ingenuity Systems, www.ingenuity.com). The software was used to identify the 
most significant molecular functions. 

"Number of functional analysis molecules related to the pathway. 

'TSfumber of functional analysis molecules related to the representative pathway subgroup. 

''Range of significance for pathway subgroups. Significance was determined Irom a single test P-value calculated using the right-tailed Fisher's exact test, 
''significance for representative pathway subgroup. 



expression profiles provides the opportunity to explore the role 
of genes with as yet unknown molecular function. This can be 
achieved by investigating interacting molecules or pathways 
the gene of interest correlates with. Our analysis identified 
several candidate genes that can serve to generate novel hypoth- 
eses and study of the function of these proteins in EAE and MS 
pathogenesis. 

We additionally identified highly significant cw-regulated 
eQTLs that do not overlap with known EAE QTLs but could 
still denote functions important for autoimmunity. Particularly 
interesting are family members of genes that associate with 
MS including Rgs4, Cyp2rl, Tmeml84a, Zfp61 and Slc27a3. 
Rgs4 is a member of the protein family of negative regulators 
of G-protein signaling (RGS) which deactivates G protein 
coupled receptor signaling (59,60). This signaling was asso- 
ciated with gene network B that correlated negatively to EAE se- 
verity phenotypes. The RGS protein family is highly regulated at 
the transcriptional level and up-regulation of certain RGS pro- 
teins decreases immune cell migration and reduced chemokine- 
dependent calcium flux. Interestingly, these are emerging as 
potentially important drug targets. Our functional analysis 
associated Rgs4 with cell proliferation, cell migration, integrin 
and endocytosis signaling. Another family member of RGS, 
Rgsl, is associated with MS and is suggested to have a role in 



immune cell regulation. Cyp2rl associated with B cell functions 
and signaling and Tmeml84a with migration of cells (Supple- 
mentary Material, Table S3 and S4). Zfp61 is associated with 
the organization of the ER and autophagy of cells. The SLC 
superfamily of transporters is responsible for the uptake of 
amino acids, peptides, ions, hormones and drugs (61) and SLC 
mRNA levels are dysregulated in inflammatory bowel disease 
patients (62). Further studies on the functions and pathways 
affected by MS family member genes identified in this study 
could add to the understanding of the genetic contribution of sus- 
ceptibility and/or pathogenesis of autoimmune diseases. 

Several c/5-eQTLs are genes with a recognized role in the 
immune system. Ccr6 has previously been demonstrated to 
play an important role in the mechanism of autoreactive lympho- 
cyte priming and migration to the efferent lymphatics (63), 
implicated in the migration of lymphocytes to the CNS (64) 
and has been associated with Crohn's disease and rheumatoid 
arthritis (47,65). Tgfb2 is a potent immunosuppressive cytokine 
that ameliorates EAE and has been associated with remissions in 
MS. However, cases of nephrotoxicity in clinical trials (66) ruled 
it out as a treatment alternative for MS patients. Ifitm6 is part of a 
family of genes that are induced by IFNs (67) and encode cell 
surface proteins that modulate cell- cell adhesion and cell differ- 
entiation. 
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Figure 3. Gene network associated with GPCR signaling correlates negatively with EAE. Expression traits with neighboring locations are drawn next to each other. 
eQTLs are indicated by their rat chromosome location on a second line and cu-eQTLs are indicated in white. eQTLs within EAE QTLs are indicated by a grey box with 
the QTL name listed in parentheses. Predicted hub genes are presented in orange and the black arrow heads indicate a directed interaction. Plus and minus indicate 
association to EAE phenotype. Genes with a higher disease association are placed to the right, whereas causative genes are placed toward the left. 



Genes can exert a control over the expression of another gene 
or a set of genes, which are referred to as frawj-eQTLs. In our 
study, the sensitivity to detect fra«5-eQTLs has greatly been 
increased by the use of a BC population originating from two 
inbred strains and Affymetrix Gene arrays that provide compre- 
hensive genome coverage with multiple probes per each tran- 
script. More than 70% of eQTLs displayed regulation in trans 
at a genome-wide corrected P-value of <0.05, which conforms 
to previous studies in inbred rodent strains (14, 68). The 
?raw5-band on RN017 overlaps Eae23a, which we previously 
demonstrated to control EAE severity in a congenic strain 
(26). Although this trans-band was not significantly enriched 
for functional pathways, we observed regulation of many 
small RNA genes such as nucleolar and spliceosomal RNAs, 
previously implicated in epigenetic regulation (31) and alterna- 
tive splicing (32), respectively. TVan^-regulation can denote a 
genomic location of a master regulator of transcription (16). A 
transcription factor, Zebl (69), known to regulate the IL-2 
pathway (70, 71) that is crucial in MS and EAE, is encoded in 
Eae23a. In addition, we have previously demonstrated that the 
balance between splice variants of Zebl driven by Eae23a 
could influence the regulation of EAE (26). 

With the recorded clinical and molecular parameters during 
the full EAE experiment (susceptibility, severity and molecular 
phenotypes), we defined clusters of gene networks that had a 
direct relation to clinical characteristics of disease. Thus, al- 
though our study in chronic stage of EAE does not necessarily 
address early molecular events that control susceptibility, a pos- 
sibility to correlate gene expression with clinical phenotypes 
provides a major benefit in defining genes and pathways that 
are related to the disease (and not just a consequence of 



genetically driven expression differences between the strains). 
Moreover, defining pathways that regulate chronic stage of 
disease might have better translational prospects considering 
that all interventions in MS occur long after the onset of 
disease. The identified genetically driven pathways that regulate 
susceptibility to autoimmune disease primarily involved T cell 
activation pathways, G protein coupled receptor and IFN signal- 
ing. The most positively correlated network to EAE susceptibil- 
ity and severity, gene network A (Fig. 2), associated with 
molecular fianctions including T cell-mediated immune mechan- 
isms. EAE and MS are described as being T cell-mediated since 
adoptive transfer of T cells can induce disease, T cell infiltrates 
are evident in EAE and MS lesions (72-75) and that the major 
genetic determinant of both MS and EAE is the HLA/MHC. 
Thus our unbiased way of defining functions central in EAE 
conform to previous knowledge about central mechanisms in 
autoimmune inflammation, indicating that these mechanisms 
are partially genetically regulated, as implicated in the MS 
GWAS. 

Themis, Lefl, Lck, Satbl, Crtam and Itk were identified in 
network A, which strongly correlated with EAE susceptibility 
and severity and associated with molecular functions including 
T cell-mediated immune mechanisms. C/s-regulated Themis 
act through TCR signaling and has recently been reported to 
control Treg functions and susceptibility to intestinal inflamma- 
tion (76). In network A, Themis was determined to directionally 
interact with the lymphoid enhancer-binding factor 1 {Lefl) that 
binds to a functionally important site in the TCR-a enhancer. In 
turn, Lefl was in our co-expression analysis predicted to interact 
with Lck, known to play a key role in the TCR signal transduction 
pathway. Lefl also interacted with Satbl, a genome organizer 
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Figure 4. Gene network associated with IFN patliway signaling. Expression traits with neigliboring locations are drawn next to each other. eQTLs are indicated by their 
rat chromosome location on a second line and cu-eQTLs are indicated in white. eQTLs within EAE QTLs are indicated by a grey box with the QTL name listed in 
parentheses. Predicted hub genes are presented in orange and the black arrow heads indicate a directed interaction. Plus and minus indicate association to EAE pheno- 
type. Genes with a higher disease association are placed to the right, whereas causative genes are placed toward the left. 



that regulates chromatin structure and gene expression, known to 
be crucial for the phenotype and function of Tregs (77). 

The T-cell-specific HMG box transcription factor family 
member Tcf7, also part of network A, is a transcriptional activa- 
tor for genes involved in immune regulation (78). The encoded 
protein can bind an enhancer element and activate the CD3E 
gene that is part of the T cell receptor- CD3 complex. We also 
found a directed interaction between cytotoxic and a regulatory 
T cell molecule { Crtam), found to influence the adaptive immune 
response (79), and the IL-2-inducible T cell kinase (Itk), which is 
thought to play a role in T cell proliferation and differentiation 
(80). Lefl and the hub gene phosphodiesterase 3B (PdeSb) inter- 
acted with Cd40l expressed on activated T cells. Cd40l is a cost- 
imulatory molecule and induces activation in APC by binding to 
Cd40. 

Utilizing information about clinical EAE development 
(disease status, severity and onset) enabled predictions regarding 
the direction of interactions demonstrating a causal relationship 
associated with the development of a clinical phenotype. 



Additionally, our analyses can serve to reveal novel players 
that influence specific molecular pathways. For instance, in the 
gene network associated with T cell functions, the co-expression 
analysis predicted that Mfsd4 would interact, in a directed 
causal-effect relationship PdeSb. Pde3b is an enzyme that 
hydrolyzes cAMP for cell metabolism and has been shown to 
be regulated downstream of Foxp3, by direct binding, to 
support Treg homeostasis and lineage stability (81,82). 

We identified a gene network highly associated with IFN sig- 
naling pathways which included STAT transcription factors, 
interferon regulatory factors and IFN-induced downstream sig- 
naling molecules. A first-line therapy for MS is IFN-|3, a type I 
IFN, although the exact molecular mechanism of action of the 
drug is still poorly understood. Mxl and Mx2, part of the 
network, can be used as biomarkers for monitoring response to 
interferon therapy in patients (83). Mx2 expression is regulated 
in cis on RNOll, and in IPA we found Mx2 to associate with 
IFN signaling and activation of IRF by cytosolic pattern recogni- 
tion receptors, attraction of NK cells, homing of T lymphocytes 
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and cell death of phagocytes (Supplementary Material, Tables S3 
and S4). In gene network C, we also recorded the nuclear body 
proteins (Pml and SpllO) that regulate genes involved in the 
immune response (84). Pml acts as a transcriptional co-activator 
with p53 (85). p53 has previously been reported to interact with 
Wwox (86), a cw-eQTL on RN019. WW domain-containing pro- 
teins are involved in protein degradation, transcription and RNA 
splicing. The network associated with IFN signaling included 
several T cell chemoattractants {Xcll, Cxcl9 and CxcllO). Xcll 
is a cytokine that attracts T cells and is expressed in activated 
thymic and peripheral blood CD8''' T cells (87, 88). Cxcl9 and 
CxcllO are induced by IFN-7 and execute their chemotactic func- 
tions by interacting with the chemokine receptor CxcrS. 

In conclusion, our well powered and unique study of the 
non-MHC iniluences on gene expression generates hypothesis 
of molecular pathways that are genetically regulated and predis- 
pose for EAE. We detected cw-regulated transcripts, both in 
EAE QTLs and outside, which may play important roles in regu- 
lating key mechanisms in both EAE and human disease owing to 
a large overlap of genes and functions. These can serve to gener- 
ate novel hypotheses useful in further dissecting pathogenic mo- 
lecular mechanisms that are dysregulated during chronic 
autoimmune inflammation. In addition, the gene network ana- 
lyses denoted genes linked to functional pathways and described 
their importance in regulating clinical traits of disease. We dem- 
onstrate how MS risk genes translate well to EAE and efforts 
should now be put into combining clinical and experimental re- 
search to evaluate candidate genes and to functionally study their 
mechanisms of action. Ultimately, these integrated findings can 
provide insight to possible diagnostic and prognostic biomarkers 
or potential therapeutic interventions for autoimmunity. 



MATERIALS AND METHODS 

Experimental subjects 

DA rats were originally obtained from the Zentralinstitut fiir 
Versuchstierzucht (Hannover, Germany) and MHC-identical 
PVG.lavl from Harlan UK Limited (Blackthorn, UK). Colonies 
have thereafter been established at Karolinska Institutet (DA/ 
Kini and PVG.lavl/Kini). A (DA x PVG.lavl) x DA BC 
population was established. In short, to create the Fl generation, 
four breeding pairs with DA female founders were established. 
The N2 generation was bred from eight breeding pairs, with 
DA females or males crossed to Fl males and females, respect- 
ively. FourN2 litters, 421 BC rats (213 females and 208 males), 
were subjected to MOG-EAE in four separate experiments, re- 
ferred to as experimental sets. Out of 42 1 immunized BC rats, 
splenic RNA was extracted from 347 rats (182 females and 
165 males) that survived until day 35 p.i. All 421 BC rats were 
used to map EAE pQTLs, which are indicated in Supplementary 
Material, Table S5 and the details of mapping will be reported in 
a separate study (Stridh et al., unpublished data). To minimize 
variability introduced by experimental sets and increase 
mapping power, we decided to perform genome-wide expres- 
sion analysis in one sex; male rats were selected due to less vari- 
ability between experimental sets. Out of 165 males, 150 
selected rats represent similar susceptibility and severity, and 
all breeding pairs and experimental sets as the full set of 165 
rats. All animals were bred and housed at the Karolinska 



University Hospital (Stockholm, Sweden) in 12h light/dark- 
and temperature-regulated rooms in polystyrene cages contain- 
ing aspen wood shavings where they had free access to standard 
rodent chow and water. Rats were tested according to a health- 
monitoring program at the National Veterinary Institute (Statens 
Veterinarmedicinska Anstalt, SVA) in Uppsala, Sweden. 



Induction and clinical evaluation of EAE 

Recombinant MOG, amino acids 1-125 from the N-terminus, 
was expressed in Escherichia coli and purified to homogeneity 
by chelate chromatography as previously described (89). 421 
BC rats, between 10-14 weeks of age, were anesthetized with 
isofluorane (Forene, Abbott Laboratories, Abbot Park, IL, 
USA) and immunized with a single subcutaneous injection at 
the dorsal tail base with 200 |jl1 of inoculum containing rMOG 
(12.5-15 |JLg/rat in females and 25-30 |JLg/rat in males titrated 
in order to achieve similar disease severity/induction) in saline 
emulsified in a 1:1 ratio with incomplete Freund's adjuvant 
(Sigma Aldrich, St Louis, MO, USA). Rats were monitored 
daily for weight and clinical signs from day 8 until 35 p.i. The 
clinical score was graded as follows: 0, no clinical signs of 
EAE; 1, tail weakness or tail paralysis; 2, hind leg paraparesis 
or hemiparesis; 3, hind leg paralysis or hemiparalysis; 4, tetra- 
plegia or moribund; 5, death. Clinical parameters were assessed 
and used in the analysis: EAE, incidence of EAE; ONS, onset of 
EAE (day p.i. of first clinical sign); DUR, duration of EAE 
(number of days subjects showed clinical signs of EAE); 
MAX, maximum EAE score during the experiment; SUM, 
sum of all scores during the experiment; s35, EAE score on 
day 35 post immunization (p.i.) (last day of experiment when 
spleen tissue was collected for expression analysis); WL, 
weight loss (calculating the percentage of weight loss between 
lowest weight throughout the experiment and day 8 p.i.). At 
day 35 p.i., 347 rats were sacrificed and spleens were collected, 
snap-frozen and stored in — 70°C until use. Spleens from the 
remaining 74 animals were not collected due to decease prior 
to the end of experiment. In the BC experiment, 60.4% of all 
rats were affected with EAE (average maximum clinical 
score = 2.46). Of individuals still included in the experiment 
at day 35 p.i., 6 1.8% were still affected at the end of experiment 
(average clinical score at day 35 p.i. score = 2.05). Anti-MOG 
antibodies were measured at the onset (day 12 p.i.) and chronic 
phase of EAE (day 35 p.i.) as previously described (90). 



DNA isolation and genotyping 

Genomic DNA was prepared from tail tips of the BC population 
as described previously. Information about polymorphic micro- 
satellite markers for the BC (118 markers evenly spaced 
throughout the genome with an average inter-marker distance 
of 20 cM) was retrieved from Ensembl Genome Database (http:// 
www.ensembl.org v. 50-62) (91). DNA amplification was per- 
formed with PCR using forward primers end-labeled with a 
fluorescent dye (VIC, NED, FAM or PET) with products run 
on ABI 3730 capillary sequencer and analyzed with GeneMap- 
per v3.7 (Applied Biosystems, Foster City, CA, USA). Primers 
were obtained from Eurofins MWG Operon (Ebersberg, 
Germany), Applied Biosystems or Proligo (Paris, France). All 
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genotypes were evaluated manually and quality assessed by two 
independent observers. 

Expression analysis 

One hundred milligrams of splenic tissue from each of 150 BC 
rats was disrupted using Lysing Matrix D tubes (MP Biomedi- 
cals, Irvine, CA, USA) in a FastPrep homogenizer (MP Biome- 
dicals). mRNA was extracted using an RNeasy mini kit 
according to the manufacturers protocol (Qiagen, Hilden, 
Germany), including on-column DNA digestion. RNA concen- 
tration and purity was determined through measurement of 
A260/A280 ratios with a NanoDrop ND-1000 Spectrophotom- 
eter (NanoDrop Technologies, Wilmington, DE, USA). For 
150 BC male rats selected to represent all breeding pairs and 
sets, whole-genome expression profiles were determined with 
Affymetrix Rat Gene 1.0 ST Array (Affymetrix, Santa Clara, 
CA, USA). Confirmation of RNA quality was assessed using 
the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa 
Clara, CA, USA). Target labeling, array hybridization and 
washing and staining were performed as described in the Gene- 
Chip Whole Transcript (WT) Sense Target Labeling manual 
(http://www.affymetrix.com). Arrays were scanned in the Gen- 
eChip Scanner 3000 7G (Affymetrix). Affymetrix cell intensity 
(CEL) files from scanning were analyzed with Affymetrix Ex- 
pression Console (EC, version 1.1). The following settings 
were used: Summarization: probe logarithmic intensity error 
(PLIER) described in the Guide to probe logarithmic intensity 
error (PLIER) estimation (http://www.affymetrix.com); Back- 
ground correction: PM-GCBG and Normalization: Global 
Median. The microarray data are available in MIAME- 
compliant (minimal information about a microarray experi- 
ments) format at the ArrayExpress Database (http://www.ebi. 
ac.uk/arrayexpress) (92) under accession code E-MTAB-784. 

Linkage and statistical analysis 

For the BC, genome-wide linkage analysis was performed for 
27342 expression traits with a set of 1 1 8 genetic markers using 
the Haley-Knott regression model in R/qtl version 1.14 (93), 
to generate maximum base 1 0 logarithm of the likelihood ratio 
of odds (LOD) scores. Logarithmic values of all molecular quan- 
titative traits (transcript expression) were used as an input for 
R/qtl and analyses were performed independently for several 
clinical phenotypes to act as covariates as well as single QTL 
scans for all covariates and pairs of covariates. Genome-wide 
significance (P- value <0.05) was generated with 1000 permuta- 
tion tests (94). Physical locations of probe sets were obtained 
from Ensembl version 55 or from Affymetrix (http://www.a 
ffymetrix.com) and the transformation of cM to Mb was deter- 
mined with R/qtl. The physical positions (Mb) were retrieved 
from Ensembl, either directly or by sequence similarity searches 
with the oligo sequences. For arbitrary cM, the density of the 
markers was considered sufficient (95) to allow a linear interpol- 
ation between the physical positions of flanking markers. Com- 
putations and their interpretation were orchestrated with the 
interactive QTL system (TIQS. http://tiqs.it), where data are 
available. eQTLs were defined as cis when there was no 
genetic marker between the peak of the linkage score and the 
chromosomal region coding the transcript. Other eQTLs were 



considered being regulated in trans. Genomic locations were 
obtained from genomic assembly Rnor3 .4. For pathway analysis 
in IPA, significance was determined with the right-tailed 
Fisher's exact test and adjusted significance using the Benja- 
mini-Hochberg correction. Venn diagrams are based on the 
gplots R package version 2.11.0.1 (http://CRAN.R-project. 
org/package=gplots). 

Due to the polygenic nature of EAE, we used a multiple-QTL 
model, i.e. forward selection followed by backward elimination 
in R/qtl version 1.14(93), to identify EAE QTLsin421 BCrats, 
encompassing 150 rats used for eQTL analysis (Supplementary 
Material, Table S5). Similar results were obtained using Hailey- 
Knott regression (data not shown). EAE QTL analysis was also 
performed in 150 rats used for eQTL analysis using the Haley- 
Knott regression method in R/qtl version 2.12 (93). Nominal 
P- values were used to report EAE QTLs in 150 rats due to 
prior evidence from 421 BC rats and additional independent 
studies (Supplementary Material, Table S5). 

Functional association analysis and network analysis 

Molecular functions and canonical pathways were evaluated 
with the ingenuity pathways application (Ingenuity Systems, 
www.ingenuity.com). The functional and canonical pathways 
(disease-specific pathways were not included) analysis identi- 
fied the molecular functions and pathways, respectively, from 
the Ingenuity Knowledge Base that was most significant to the 
uploaded data set. For analysis of transcripts correlated to 
cw-eQTLs (Pearson's correlation coefficient r > 0.40 or 
r < — 0.40), the experimental set 1 was excluded due to signifi- 
cant deviation from other sets (the induction dose of 30 |JLg/rat 
was used for set 1 males, whereas other sets were induced with 
25 fJLg/rat). 

Gene networks were provided by TIQS, created on the basis of 
WGCNA (35). The Pearson correlation coefficient was used to 
determine the association of a gene network with clinical EAE 
phenotypes. Gene network analysis in combination with the clin- 
ical parameters measured for each individual in the EAE experi- 
ment was used to predict a causal relationship between genes. 
Directions indicate the gene that is observable to differ in its ex- 
pression level from healthy individuals after longer disease dur- 
ation. Directed associations were determined with the use of 
clinical parameters for all predicted gene interactions and hub 
genes in the eQTL data as inferred from an association of the 
gene expression levels with the onset of the disease phenotypes 
(Gupta etal. , unpublished data). The arrows point to the gene that 
changes its expression at a later stage of the disease. To deter- 
mine the direction, the expression levels were compared with 
the onset and the severity. From the onset, the duration of the 
disease is derived. Every gene is assigned an average score 
across all individuals to reflect if the expression levels are chan- 
ging with the duration of the disease, with the strength of the 
phenotype contributing to the scoring. The directions follow 
the gradient of that score. 
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